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Abstract 

We study the effect of coadsorption of CO and O on a Ziff-Gulari-Barshad (ZGB) 
model with CO desorption (ZGB-d) for the reaction CO -|- O — >• CO2 on a catalytic 
surface. Coadsorption of CO at a surface site already occupied by an O is introduced by 
an Eley-Rideal-type mechanism that occurs with probability p, < p < 1. We find that, 
besides the well-known effect of eliminating the second-order phase transition between the 
reactive state and an 0-poisoned state, the coadsorption step has a profound effect on the 
transition between the reactive state and the CO-poisoned state. The coexistence curve 
between these two states terminates at a critical value kc of the desorption rate k which 
now depends on p. Our Monte Carlo simulations and finite-size scaling analysis indicate 
that kc decreases with increasing values of p. For p = 1, there appears to be a sharp phase 
transition between the two states only for k at (or near) zero. 



PACS numbers: 82.65. +r, 68.43.De, 05.50.+q, 05.10.Ln 



I. INTRODUCTION 



Catalytic reactions on surfaces are very complex, and their comprehensive un- 
derstanding is vital, both for its intrinsic scientific interest and for its immense 
technological applications.-i^ The scientific interest is due to the emergence of a rich 
and complex variety of phenomena, including chaotic behavior, bistability, critical 
phenomena, out-of-equilibrium phase transitions, etc.i'^ Due to its crucial role in 
industrial applications, the oxidation of CO on a transition-metal catalyst is one of 
the most studied reactions. In a pioneering work, Ziff, Gulari and Barshad (ZGB)^ 
proposed a deceptively simple model to describe some kinetic aspects of the reac- 
tion CO -|- O ^ CO2 on a catalytic surface in terms of a single parameter which 
represents the probability that the next molecule arriving at the surface is a CO, 
i.e.,., it is proportional to the partial pressure of CO. The model exhibits two kinetic 
phase transitions, one continuous, at y = yi, between an oxygen-poisoned phase and 
a reactive one, and a discontinuous one, at y = y2, between a reactive phase and a 
CO poisoned one. Unfortunately, there are important aspects of this catalytic reac- 
tion that are not reproduced by this simple model. Transitions between states of low 
and high CO coverage, 6*00 (the fraction of surface sites occupied by CO), have been 
observed experimentally- At low temperatures, as y increases, there is a discontinu- 
ous drop in the CO2 production rate; while above a certain critical temperature the 
discontinuity disappears, and the CO2 production decreases continuously. This type 
of behavior can be reproduced by including a CO desorption rate, k, the so-called 
ZGB-d model.-i^i^ Another feature that is erroneously predicted by the original ZGB 
model is the continuous phase transition at ?/i > 0, since experiments show that the 
reaction rate increases as soon as the CO concentration departs from zero. - The real 
system does not present an oxygen-poisoned state because oxygen does not impede 
the adsorption of C0.-»^ Several authors have shown that this transition can be 
eliminated by adding to the Langmuir-Hinshelwood (LH) mechanism that defines 
the original ZGB model an Eley-Rideal (ER) step that allows a reaction between 
CO molecules in the gas phase and adsorbed O atoms on the surface.— However, 
there are no experimental data indicating that such a reaction between free CO and 
adsorbed oxygen occurs, but there is experimental evidence that suggests that under 
certain conditions, it is possible that a CO molecule can be weakly adsorbed at a 
site already occupied by an oxygen atom.J^ Experiments based on photoemission 
electron microscopy have shown that CO can be adsorbed on a saturated O sur- 
face.— This coadsorption is a prerequisite for the reaction to happen, and it can be 
simulated by an ER-type mechanism that occurs with a certain probability. Some 
previous work on the ZGB model with coadsorption indicates that the production 
rate of CO2 increases linearly with y (in the region of low y), and that the continuous 
phase transition from the reactive state to the O poisoned state disappears.— 'i^'^^'i^ 
In these works it was assumed that the addition of the ER-type step does not alter 
the nature of the transition at 2/2 ? but simply shifts it to a lower value of 7/2- 

We believe that the effect of the ER-type step on the ZGB-model deserves a 
more detailed analysis. In this work we therefore explore what happens when a 
coadsorption mechanism is added to the ZGB model with CO desorption. For the 
ZGB-d model, there is a distinction between the high and low CO-coverage phases 
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only for desorption rates k below some critical value kc, while above k^ the CO 
coverage varies smoothly with y. Thus, the transition value 1/2 becomes a function 
of k, corresponding to a coexistence curve 1/2 (fc) that terminates at the critical point 
1/2 (fee) •—'^'^ The value of kc has been estimated by several authors. An earlier 
work based on fractal scaling of the interface between phases indicates that the 
first-order phase transition disappears for k > kc ^ 0.039.^ More recent estimates 
based on the probability distribution for the CO coverage (histograms), the fourth- 
order cumulant, and finite-size scaling theory give the values kc = 0.040,— and 
kc = 0.039," respectively. In this work we are particularly interested in studying 
how the nature of the transition at 1/2 changes when the ER-type mechanism is 
added with probability p to the ZGB-d model. To the best of our knowledge, such 
an analysis has not been performed previously. Our study is based on kinetic Monte 
Carlo simulations and finite-size scaling analysis of the data. 

The rest of this paper is organized as follows. In Sec. II, we define the model 
and describe the Monte Carlo simulation techniques used. In Sec. Ill A, we present 
and discuss the numerical results obtained for the case of k > 0. In Sec. Ill B, 
we perform a finite-size scaling analysis for the particular case of the standard ZGB 
model with an Eley-Rideal step {k = 0,p = 1). Our conclusions are summarized in 
Sec. IV. 

II. MODEL AND SIMULATION 

The ZGB model with desorption and coadsorption is simulated on a square lattice 
of linear size L that represents the catalytic surface. A Monte Carlo simulation 
generates a sequence of trials: CO or O2 adsorption with probability 1 — k and CO 
desorption with probability k. For the adsorption a CO or O2 molecule is selected 
with probability y and 1 — y respectively.^'^*^ These probabilities are the relative 
impingement rates of both molecules and are proportional to their partial pressures. 
The algorithm works in the following way. A site i is selected at random. For 
desorption, if i is occupied by CO the site is vacated; if not, the trial ends. For 
adsorption, if a CO molecule is selected it can be adsorbed at an empty site i if none 
of its nearest neighbors are occupied by an O atom. Otherwise, one of the neighbors 
occupied by O is selected at random and removed from the surface, leaving i and the 
selected neighbor vacant. This move simulates the CO + O ^ CO2 surface reaction 
following the adsorption of CO. If the selected site is occupied by an O atom the CO 
is coadsorbed with probability p, and a CO2 molecule is liberated leaving behind an 
empty site on the surface. This represents the ER reaction. O2 molecules can be 
adsorbed only if a pair of nearest-neighbor sites are vacant. If the adsorbed molecule 
is selected to be O2, a nearest neighbor of i, j, is selected at random, and if it is 
occupied the trial ends. If both i and j are empty, the trial proceeds, and the O2 
molecule is adsorbed and dissociates into two adsorbed O atoms. If none of the 
remaining neighbors of i is occupied by a CO molecule, the one O atom is located 
at i, and if none of the neighbors of j is occupied by a CO molecule, then the other 
O is located at j. If any neighbors of i are occupied by a CO, then one is selected at 
random to react with the O at i such that both sites are vacated. The same reaction 
happens at site j if any of its neighbors are filled with a CO molecule. This process 
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mimics the CO + O ^ CO2 surface reaction following O2 adsorption. A schematic 
representation of this algorithm is given by 



CO(g) + S ^ CO (a) 

02((7) + 2S ^ 20(a) 

CO(a) + 0(a) ^ C02(g) + 2S 

CO (a) ^ CO(g) + S 

CO(g) + 0(a) ^ 0(a) + CO(c) 

CO(c) + 0(a) ^ C02(g) + S 

Here S represents an empty site on the surface, g means gas phase, a means adsorbed, 
and c means coadsorbed. The first three steps correspond to the LH mechanism, the 
fourth step to the CO desorption, and the last two steps to the ER-type mechanism, 
which occurs with probability < p < 1. If p = 0, we recover the ZGB-d model, 
if p = 1, we have the ZGB-d model with an ER step. It is important to realize 
that the ER reactions are essentially instantaneous, so that the CO cannot remain 
coadsorbed without reacting. 

For our simulations we assume periodic boundary conditions. The time unit is 
one Monte Carlo Step per Site, MCSS, in which each site is visited once, on average. 
Averages are taken over 10^ independent simulation runs for each set of parameters. 



III. RESULTS 

Starting from an empty lattice, we wait until the system reaches a steady state 
at constant partial pressure y. We calculate the fraction of sites occupied by CO 
molecules: the CO coverage, 9co] the O coverage, 9o, and the fraction of empty 
sites, 6*3- -Rco2 is the rate of production of CO2. 



A. Results for k > 

In Fig. [H we show the dependence of the coverages and the CO2 production 
rate on the CO pressure, y, for the particular case of k = 0.02 and p = 0.5. It is 
immediately seen that the inclusion of the ER-type step {p > 0) has a dual effect: one 
is to eliminate the oxygen-saturated phase and therefore the continuous transition 
at yi in the ZGB model, and the other is that the transition to a CO-saturated 
surface occurs at a lower value of the CO pressure than for p = 0. These two effects 
can be seen in more detail in Fig. [21 where we plot the CO coverage vs y for a fixed 
value of k and different values of p. As p increases, the second transition is shifted 
toward lower values of 1/2, and the CO coverage starts to increase as soon as y > 0, 
except for the case p = (the ZGB-d model), in which it starts at yi ~ 0.388. A 
closer scrutiny of Fig. [2] also reveals that the discontinuity in 6^00 at y2{k = 0.03, p) 
appears to become smaller as p increases, giving way to a continuous change for p 
below approximately 0.6. 
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Previous studies show that when CO desorption {k > 0) is added to the original 
ZGB model, the ZGB-d model, the discontinuous transition ends at a critical value 
kc , such that in this case 1/2 = 2/2(fcc)-~'^'^- To the best of our knowledge, there 
have been no studies of how the nature of the transition changes when an ER-type 
step is added to the ZGB-d model, and this is the main topic we explore in this 
paper. 

An efficient way to locate and classify phase transitions is to look at the fourth- 
order reduced cumulant of the order parameter. For ^co it takes the form,— i2ii2^i2^ 



Ul = 1- 



/i4 



where 



/in = ((^co - (^co))") = / (^co - {eco)rP{0co)d9co (2) 

Jo 

is the nth central moment of the CO coverage, and P{6co) is the probability dis- 
tribution for ^co- This probability distribution is obtained by counting the number 
of times, Ni, that the coverage falls in the intervals [0, A), [A, 2A), [1 — A, 1] 
(A = 0.01), such that ^jiVj = is the total number of MCSS. Then, the proba- 
bility. Pi, that ^co has a value in the interval [{i — l)A,iA) is 



such that 



[' P{9co)d9co-y^PrA = l. (4) 

"'0 



In Fig. [3] we show P{6co) for fixed values of k and p, and values of y in different 
regions: below and above 1/2, where the histograms are unimodal, indicating that 
the system consists of one single phase, and at 1/2 where the histograms are bimodal, 
signaling the existence of two distinct phases. At the coexistence point y2{k,p), the 
areas under both peaks are equal.— A finite-size scaling analysis would have to 
be performed to ascertain if the transition is continuous or discontinuous for any 
particular values of k and p. Such a study is planned for future work. 

The equal-area bimodal distribution corresponding to coexistence yields a posi- 
tive maximum for the cumulant ul y, fianked on either side by negative minima 
and approaching zero far away from the transition. The maxima of ul define the 
L-dependent coexistence line, y2{k,L,p). In Fig. IHwe show the dependence of ul 
on y for several values of of k and p = 0.5 (Fig. Ufa)), and p = 1.0 (Fig. Hl^b)). In the 
proximity of a first-order phase transition, the maximum value of the fourth-order 
cumulant (when L —>■ 00) takes the value of 2/3,— while for a continuous transition 
the maximum takes a lower value. For the Ising universality class, the maximum 
value of Ul at the critical point is approximately 0.61.^^ Our results strongly suggest 
that when the ER-type step is added (p > 0) the coexistence curve ends at some 
value of kc that decreases when p increases, i.e., now y2{kc,p). Figure |l](a) indicates 
that when p = 0.5, the first-order phase transition seems to terminate at some value 
of k below 0.02, i.e., 0.01 < kc{p = 0.5) < 0.02, while Fig. H^b) indicates that even 
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for the well-known ZGB model [k = 0, no desorption) with p = 1,^ it is not clear 
that the transition remains first-order. In the following we analyze this particular 
case in further detail. 

B. Results for k = 0,p = 1 

The standard ZGB model with the addition of an ER step has been analyzed by 
several authors,— '^'^'^'^ but as far as we know the nature of the transition at y2 
has not been studied in detail, and it has been assumed that, besides destroying the 
continuous transition at yi, the effect of the extra step is simply to shift the first- 
order transition to a lower value of ?/2- Our analysis of the fourth-order cumulant 
suggests that is depressed by increasing p, and in this section we perform finite- 
size scaling analysis to further investigate this point. 

The fluctuations of Oqq can be calculated in the standard way as 

XL = L\{elo)L-{eco)l), (5) 

where 

mo)L= f 0ZoPidco)deco. (6) 
Jo 

The peak positions of xl approach the infinite-system transition point with increas- 
ing L. In Fig. we show xl as a function of y for k = and p = 1 for several 
values of L. As expected, the peaks of the curves shift and increase in height with 
increasing L. It is well known that at a first-order equilibrium phase transition the 
order parameter fluctuations scale as xl ~ L'^, where d is the spatial dimension of 
the system (here d = 2)M^^ Previously we have successfully used the same scaling 
relation to identify the nonequilihrium first-order transition in the ZGB-d model.— 
In Fig. El^a) we plot \n{xY,^^) ln(L). The linear fit indicates a power-law scaling, 
such that Xl^"" ~ L'^ > where d' = 1.81 ± 0.02. We also measure the power-law expo- 
nent by a different method, which has some advantage in eliminating the effects of 
a nonsingular background term (as in xl = f + gL'^ with / and g constants). We 
write 

/lnb = d' + 0{l/lnb) (7) 

with L fixed at a relatively small value (here, L = 20), and 6 > 1. For large L 
and b, the correction term is proportional to f/{glnb), so that the exponent can be 
estimated by plotting the left-hand-side of Eq. (I7l) vs 1/ In 6 and extrapolating to 
1/ In 6 = 0, as in Fig. EJ^b). The resulting estimate is d' = 1.83 ± 0.02. These results 
are relatively far from the expected value d' = 2 for a first-order phase transition, 
in clear contrast with the exponent calculated by the same techniques for p = 0, 
in which case we found that d' ^ 2 for the transition at y2{k < kc,p = 0)M This 
analysis provides further support to our previous suggestion, based on the fourth- 
order cumulant, that the ZGB model with an ER step {p = 1) does not have a phase 
transition but only a smooth crossover for any non-zero value of k. 
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IV. CONCLUSIONS 



In this paper we have presented a study of the ZGB-d model in which an ER- 
type step is added (0 < p < 1). The extra step can be thought of as a way to 
include a surface interaction between CO and O coadsorbed at the same site. We 
found that, besides eliminating the unphysical oxygen-poisoned state and therefore 
the second-order transition of the original ZGB model, the inclusion of the ER-type 
step has a more profound effect than the mere shifting of the first-order transition at 
?/2 toward lower values, which was previously reported in the literature.— li^iii^ The 
addition of the ER-type step changes the nature of the transition between the CO2 
producing phase and the non-productive CO-poisoned phase, in such a way that 
now the coexistence curves terminate at a value of kc that decreases with p. Even 
for the well-studied case of the ZGB with an ER mechanism {p = 1), our results 
indicate that it is not legitimate to speak of a phase transition at ?/2 for any non-zero 
k. These results are summarized in Fig. [71 where we show a rough estimate of the 
coexistence curves indicating for some values of p the upper value of k for which 
there could still exist a first-order phase transition. We hope our results will be of 
use for the analysis of experiments on heterogeneous CO oxidation. 
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FIG. 2: CO coverage vs y for desorption rate k = 0.03 and several values of p. As p 
increases, 2/2 decreases. In the inset it can be seen that, except in the particular case 
p = 0, m which the system has a phase transition at yi ~ 0.388, the CO coverage starts 
increasing as soon as ?y > 0. L = 60. 




FIG. 3: Order-parameter probability distribution P(^co) for fixed k = 0.02 and p = 0.5, 
below (y = 0.5125, circles), and above (y = 0.5175, diamonds) the transition at y2, where 
the distributions are unimodal and remain in the low and high coverage phases respectively. 
At the transition point, 2/2 = 0.02, p = 0.5) ^ 0.5161 (squares), the histogram is bimodal 
and the areas under the peaks are equal, suggesting the coexistence of two different phases. 
L = 60. 
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FIG. 4: Fourth-order cumulant ul vs y for (a) p = 0.5 and (b) p = 1. The minima of ul 
correspond to the crossovers between unimodal and bimodal distributions. The maximum 
between these minima gives the coexistence point, y2{k,p, L). The dotted horizontal hues 
in the insets correspond to uj^ = 2/3, and the dashed horizontal Hnes to ujsing 

w 0.61. 

L = 60. For p = 0.5 (a), the maximum occurs at the same value of y (0.516), for which 
the areas under the two peaks of the order-parameter distribution in Fig. [3] are equal. 
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FIG. 5: Fluctuation of the order parameter, xl, vs y for k = 0, p = 1 for different system 
sizes L. 
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FIG. 6: Scaling behavior of Xl^^ for the ZGB model with an ER step {k = 0,p = 1). (a) 
InixT"") vs ln(L) and (b) In ix^l"" / xT"") /^^ (b) vs l/ln(b) with L=20. The straight lines 
are the best linear fits to the data and give Xl**^ ~ L'^ ■ From (a) d' = 1.81 it 0.02 and 
from (b) d' = 1.83 ± 0.02. 
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FIG. 7: Some points on the coexistence curve 7/2 shown vs k for different values of p. In 
each case the symbol with the asterisk indicates our highest estimate of kc{p), where the 
coexistence curve terminates at a critical point. For p = 1, kc appears to be at (or close 
to) zero. 
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